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ABSTRACT 

Aims. We try to determine the Galactic structure by comparing the observed and modeled velocities of OB-associations in the 3 kpc 
solar neighborhood. 

Methods. We made N-body simulations with a rotating stellar bar. The galactic disk in our model includes gas and stellar subsystems. 
The velocities of gas particles averaged over large time intervals (~ 8 bar rotation periods) are compared with the observed velocities 
of the OB-associations. 

Results. Our models reproduce the directions of the radial and azimuthal components of the observed residual velocities in the 
Perseus and Sagittarius regions and in the Local system. The mean difference between the model and observed velocities is AV = 3.3 
km s"'. The optimal value of the solar position angle di, providing the best agreement between the model and observed velocities 
is 6fc = 45 ± 5°, in good accordance with several recent estimates. The self-gravitating stellar subsystem forms a bar, an outer ring 
of subclass Ri, and slower spiral modes. Their combined gravitational perturbation leads to time-dependent morphology in the gas 
subsystem, which forms outer rings with elements of the Ri- and iJ2-morphology. The success of N-body simulations in the Local 
System is likely due to the gravity of the stellar /?i -ring, which is omitted in models with analytical bars. 
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I The consensus since the 1990s has been that the Milky 
. Way is a barred galaxy (see, e.g. iBlitz & Spergel I199U 

■ iBlitz et all [l993). The estimate for the size of the large- 
' scale bar has grown from initial Rhar ~ 2-3 kpc to 

current estimates Rh aj- = 3 - 5 kpc dHabing et al.L l2006t 

' ICabrera-Lavers et"aL i 120071 l2008h . The pos ition angle of the 
' bar is thought to be in the range 15° - 45° j Blitz etal j 11993 : 

Kuiiken, 1996; Weiner & Sellwoodl 119991; Benjam in at al. . 

2005; Englmaier & Gerhard, 2006; Cabrera-Lavers et al.L 12007 : 
. iMinchev et al.. ,2009) . The differences in the position angle es- 

■ timates may i ndicate that the innermost s tructure is actually a 
' triaxial bulge dCabrera-Lavers et" 120081) . On the other hand, 

this ambiguity may be partly caused by our unfavorable view- 
ing angle near the disk plane, which also hinders study of other 
, aspects of Galactic morphology. 

The suggested configurations for the spiral morphology of 
the Galaxy include models or sketches containing from two to 
six spiral arms (see e.g. Vallee, 2005; Vallee, 2008, and refer- 
ences therein). A case has also been suggested where a two- 
armed structure dominates in the old stellar population, whereas 
the g as and young stellar population exhibits a fou r-armed struc- 
ture (iLepine et alll200l '; Church well et all. '2009). In addition to 
spiral arms, there may be an inner ring or pseudoring surround- 
ing the bar, which manifests itself as the so-called 3-kpc arm(s) 
(iDame & ThaddeusLl2008HChurchwell etaLLl2009l) . Also, spec- 
ulations about a nuclear ring with a major axis of about 1.5 



kpc have been made (iRodriguez-Fernandez & CombesL lloosh . 
Different kinds of rings - nuclear rings, inner rings and outer 
rings - are often seen in the di sks of spiral galaxies, e specially if 
there is also a large-scale bar (lButa&Combeslll99d) . Thus, the 
presence of an outer ring in the Galaxy may also be considered 
plausible (lKalnaisl[T99ll) . 

Since the outer rings have an elliptic form, the broken 
outer rings (pseudorings) resemble two tightly wound spiral 
arms. Nevertheless their connection with the density-wave spi- 
ral arms is not very obvious because their formation does not 
need the spiral-shaped perturbation in the stellar disk. The 
main ingredient for their formatio n is a rotating ba r. Both 
test particle simulations dSchwarzL Il98ll: iBvrd et all 11994 
Baglev et al., 200 ^ with an analytical bar and N-body simula- 
tions (Rautiainen & Said 1 19991 |200tf), where the bar forms in 
the disk by instability, show that the outer rings and pseudor- 
ings are typically located in the region of the Outer Lindblad 
Resonance (OLR). Two main classes of the outer rings and pseu- 
dorings have been identified: the /?i -rings and /?j -pseudorings 
elongated perpendicular to the bar and the /?2-rings and Ri,- 
pseudorings elongated parallel to the bar. In addition, there is 
a combined morphologic a l type RfR'f that shows elem e nts of 
both classes dButal Iw, 'Buta & CrockerL [1991 iButaL [19951: 



2007h . 
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Bu ta & Combesl ll996; Buta et al. 

ISchwarzl ([1981 connected two main types of the outer 
rings with two main families of periodic orbits existing near 
the OLR of the bai^j Contopoulos & Papayannopoulos, 198(]|; 
IContopoulos & Grosboil Il989l) . The stability of orbits enables 
gas clouds to follow them for a long time period. The Ri- 
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rings are supported by xi(2)-o rbits (using the nomenclature of 
IContopoulos & GrosbolL Il989h lying inside the OLR and elon- 
gated perpendicular to the bar, while the /?2-rings are supported 
by ji(;i(l)-orbits situated a bit outside the OLR and elongated 
along t he bar. There is also another c onception of the ring for- 
mation. iRomero-Gomez et alj (l2007l) show that Lyapunov pe- 
riodic orbits around Li and L2 equilibrium points can lead to 
the formation of the spiral arms and the outer rings. They as- 
sociate the spiral arms emanating from the bar's tips with the 
unstable manifolds of Lyapunov orbits. This approach can be 
useful for explainin g of the motion of gas particles as well 
dAthanassoula et al.L [2009). 

Besides the bar the galactic disks often contain spiral arms, 
which modify the shape of the gravitational perturbation. In the 
simplest case, the pattern speeds of the bar and spiral arms are 
the same. In many studies this assumption has been used for con- 
structing the gravitational potential from near-IR observations 
(which represent the old stellar population better than the visual 
wavelengths). Several galaxies with outer rings have been mod- 
eled by this method, and findings are in good accordance with 
studies made by using analytical bars: the outer rings tend to be 
located near the OLR ( Salo etallll999l) . although in some cases 
they can be completely confined within the outer 4/1 -resonance, 
(lTreuthardtetalil2008h . 

A real galactic disk provides further complications, which 
can be studied by N-body models, where the bars and spi- 
ral arms are made of self-gravitating particles. In particular, 
there can often be one or more modes rotating m ore slowly 
than the bar (jS ellwood & Sparke, 1988; M asset & Ta gger. 1997; 
iRautiainen & Salolll99 9). Even if there is an apparent connec- 
tion between the ends of the bar and the spiral arms, it is no 
guarantee that the pattern speeds are equal - the break between 
the components may be seen only for a s hort time before the 
connection reappears (see Fig. 2 in ISellw ood & Sparke, 1988). 
Sometimes the bar mode can contain a considerable spiral pait 
that forms the observed sp iral, together with the slower modes 
dRautiainen & SaloL Il999h . The multiple modes can also intro- 
duce cyclic or semi-cycUc variations in the outer spiral morphol- 
ogy: outer rings of different types can appear and disappear tem- 
porarily (Rautiainen & Salo, 2000j). 

In lMeFnik & Rautiainen (20091) (hereafter Paper I), we con- 
sidered models with analytical bars. In this case the motion of 
gas particles is determined only by the bar We found that the 
resonance between the epicyclic motion and the orbital motion 
creates systematical noncircular motions that depend on the po- 
sition angle of a point with respect to the bar elongation and on 
the class of the outer ring. The resonance kinematics typical of 
the outer ring of subclass RiR'2 reproduces the observed veloci- 
ties in the Perseus and Sagittarius regions well. 

In paper I we also suggested that the two-component 
outer ring could be misinterpreted as a four-armed spiral. In 
some galaxies with the combined /?i/?2-morphology, the Ri- 
component can also be seen in the near infrar ed, but the 
j??-co mponent is usually prominent only in blue dBvrd et aI.L 
|1994|) . This could explain the ambiguity of the number of spi- 
ral arms in the Galaxy. N-body simulations confirm that the 
7? J -rings can be forming in the self-gravitating stellar subsys- 
tem, while the -rings us ually exist only in the gas component 
(IRautiainen & SaIoL l200(j). 

In the present paper we study the effect of multiple modes 
and their influence on the kinematics and distribution of gas par- 
ticles. We construct N-body models to study the influence of 
self-gravity in the stellar component on the kinematics of gas 
particles. We compare the model velocities of gas particles with 



the observed velocities of OB-associations in the neighborhood 
3 kpc from the Sun. 

This paper has the following structure. Observational data 
are considered in Sect. 2. Section 3 is devoted to models and 
describes the essential model parameters, the evolution of the 
stellar and gas components: formation of the bar and the inter- 
play between the bar and slower spiral modes. In Sect. 3 we also 
analyze the general features of the gas morphology. Section 4 is 
devoted to the comparison between the observed and modeled 
kinematics. Both momentary and average velocities of gas par- 
ticles are considered. The influence of the bar position angle Oi, 
on the model velocities is also investigated in Sect. 4, as are the 
evolutionary aspects of kinematics. Section 5 consists of conclu- 
sions and discussion. 



2. Observational data 

We have compared the mean residual velocities of OB- 
associations in the regions of intense star formation with those 
of gas particles in our models. These regions pr actically coincide 
with the stellar-gas complexes identified by Efremov & Sitnikl 
( 1988). The residual velocities characterize the non-circular mo- 
tions in the galactic disks. They are calculated as differences 
between the observed heliocentric velocities (corrected for the 
motion to the apex) and the velocities due to the circular rotation 
law. W e used the list of OB-associations by Blaha & Humphre^sl 
()198J, the line-of-sight velocities (Barbi er-Brossat and Figo nL 
2000), and proper motions ( Hippa rcosi 1 19971; Ivan Leeuwem 
2007) to calculate their median velocities along the galactic 
radius-vector, Vr, and in the azimuthal direction, Vg. Figure 
[1] shows the residual velocities of OB-associations in the re- 
gions of intense star formation. It also indicates the grouping of 
OB-associations into stellar-gas complexes. For each complex 
we calculated the mean residual velocities of OB-associations, 
which are listed in Table [Tj Positive radial residual velocities 
Vr are directed away from the Galactic center, and the positive 
azimuthal residual velocities Vg are in the sense of Galactic rota- 
tion. Table [1] also contains the rms errors of the mean velocities, 
the mean Galactocentric distances R of OB-associations in the 
complexes, the corresponding intervals of galactic longitudes I 
and heliocentric distanc es r, and names of OB-ass ociations the 
region includes (see also lMel'nik & DambisLl2009l) . 

The Galactic rotation curve derived from an analysis of 
the kinematics of OB-associations is nearly flat in the 3- 
kpc solar neighborhood and corresponds to th e linear veloc- 
ity at the solar distance of 0o - 220 km s ' dMel'nik et all 
I2OOIL iMel'nik & DambisL l2009l) . The nearly flat form of flie 
Galactic rotation curve was found in many other studies 
(Burton & Gordon, 1978; Clemems, 1985; 'Brand & Blitz, 1993; 
Pont. 1994: Dambis et aL. 1995.: Russeil .2003.: Bobvlev at al . 
1^0% . 

We adopted the Galactocentric distance of the Sun to be 
Rq ^ 7.5 kpc (Rastorguev et al., 1994; Dambis e taTl 119951; 
iGIushkova et aI.L 119981 and other papers), which is consistent 
with the so-calle d short distance scale for classical Cepheids 
dBerdnikovetaUl2000 ). 



3. Models 

3.1. The model parameters 

We made several N-body models, which satisfy "broad observa- 
tional constraints": the rotation curve is essentially flat and the 
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Fig. 1. (a) The residual velocities of OB-associations projected on to the galactic plane. It also shows the grouping of OB-associations 
into regions of intense star formation, (b) The mean Vr- and Vg- velocities of OB-associations in the stellar-gas complexes. The 
X-axis is directed away from the galactic center, and the Y-axis is in the direction of the galactic rotation. The Sun is at the origin. 

Table 1. Observed residual velocities of OB-associations in the stellar-gas complexes 



Region 


R 


obs 

km s" 


^0 obs 


/ 


r 


Associations 




kpc 


km s"' 


deg. 


kpc 




Sagittarius 


5.6 


+9.9 ± 2.4 


-1.0± 1.9 


8-23 


1.3-1.9 


Sgr OBI, 0B7, OB4, Ser OBI, 0B2, 
Set 0B2, 0B3; 


Carina 


6.5 


-5.8 ±3.3 


+4.7 ± 2.2 


286-315 


1.5-2.1 


Car OB 1 , 0B2, Cru OB 1 , Cen OB 1 , 

Coll 228, Tr 16, Hogg 16, NGC 3766, 5606; 


Cygnus 


6.9 


-5.0 ±2.6 


-10.4 ± 1.4 


73-78 


1.0-1.8 


Cyg OBl,OB3, 0B8, OB9; 


Local System 


7.4 


+5.3 ±2.8 


+0.6 ±2.5 


0-360 


0.1-0.6 


Per 0B2, Mon OBI, Ori OBI, Vela OB2, 
Coll 121, 140, Sco 0B2; 


Perseus 


8.4 


-6.7 ± 3.0 


-5.9 ± 1.5 


104-135 


1.8-2.8 


Per OBI, NGC 457, Cas OB8, OB7, OB6, 
OB5, 0B4, 0B2, OBI, Cep OBI; 



size of the bar is acceptable. From these models we have chosen 
our best-fitting case, which we describe here in more detail. 

The rotation curve of our best-fitting model is illustrated in 
Fig. 121 In the beginning, the rotation curve is slightly falling 
in the solar neighborhood, but the mass rearrangement in the 
disk during the bar formation makes it rise slightly. We scaled 
the simulation units to correspond to our preferred values of 
the solar distance from the Galactic center and the local cir- 
cular velocity. This also gives the scales for masses and time 
units. However, in the following discussion we will use simula- 
tion time units, one corresponding to approximately 100 million 
years, and the full length of the simulation is 6 Gyrs. 

The bulge and halo components are analytical, whereas the 
stellar disk is self-gravitating. The bulge is represented by a 
Plummer sphere, mass Mi„/ge = LITxIO^'Mq, and scale length 
Rbuige = 0.61 kpc. The dark halo was included as a component 
giving a halo rotation curve of form 



V(R) 



(1) 



where Vmax - 210 kms" is the asymptotic maximum on the 
halo contribution to the rotation curve and Rc -1.6 kpc the core 
radius. 



The N-body models are two-dimensional, and the gravita- 
tional potential due to self-gravitating particles is calculated by 
using a logarithmic polar grid (108 radial and 144 azimuthal 
cells). The N-body code we used has been written by H. Salo (for 
more details on the code, see Salo, 1991; Salo & Laurikaine^ 
l2000h . The value of the gravitation softening is about 0.2 kpc 
on the adopted length scale. The mass of the disk Mdisk = 
3.51 X 10"' Mq. 

The disk is composed of 8 million gravitating stellar parti- 
cles, whose initial distribution is an exponential disk reaching 
about 10 scale lengths. The disk and halo have nearly equal con- 
tribution to the rotation curve at the solar distance. The initial 
scale length of the disk was about 2 kpc, but after the bar for- 
mation, it forms a twin profile disk: the inner profile becomes 
steeper and the outer profile shallower, and the exponential scale 
length corresponds to about 3 kpc outside the bar region. The 
initial value of the Toomre-parameter Qt was 1.75. 

The gas disk was modeled by inelastically colliding test par- 
ticles as was done in Paper I. The initial velocity dispersion of 
the gas disk was low, about 2 km s"', but it reached typical val- 
ues in the range 5-15 km s ' during the simulation. If col- 
lisions are omitted, the velocity dispersion of the test particles 
rises much higher into the range 25 - 50 km s"' . The model used 
in the kinematical analysis contains 40 000 gas particles initially 
distributed as a uniform disk with an outer radius of 9.2 kpc. 
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Fig. 2. The rotation curve (solid line) of the N-body model at T = (left) and at T = 55 (right). The contributions from the bulge 
(dash-dotted Une), disk (dashed Une) and halo (dotted line) are also indicated. 



stars, m=2, T=5-15 Stars, m=2, T=20-30 




2 4 6 8 10 12 2 4 6 8 10 12 

R [kpc] R [kpc] 



Fig. 3. The amplitude spectra of the relative density perturbations in the model disk. The frames show the amphtude spectra of 
the stellar or gas component at various times (indicated on the frame titles). The contour levels are 0.025,0.05,0.1,0.2,0.4, and 0.8, 
calculated with respect to the azimuthal average surface density at each radius. The continuous lines show the frequencies D. and 
D. + K/m, and the dashed curves indicate the frequencies D. ± k/4 in them = 2 amplitude spectrum. 
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3.2. Evolution of the stellar component 

The inner regions quickly develop a small spiral (at T ~ 2.5), 
which then evolves to a clear bar (T ~ 5). Its original pattern 
speed Q.h is about 80 km s^'kpc"', meaning that when it forms 
it does not have an Inner Lindblad Resonance (ILR). In its early 
phase the bar slows down quite quickly (Q/, ^ 60 km s^'kpc"' 
at r = 10), but the deceleration rate soon settles down: Q/, « 
54 km s^'kpc"' at T = 20 and Q/, a 47 km s^'kpc"' at T = 55. 
In this model the bar's slowing down is accompanied by its 
growth, and the bar can always be considered dynamically fast 
(see e.g. lDebattista & Sellwo od'. 200 0J. Using the same method 
to determine the bar l ength a s Rautiai nen et al.l (2008) (a modifi- 
cation of one used bv' Erwinl 12005 ). we get Ri,a,- = 4.0 + 0.6 kpc 
at T=55 and Rcu/Rhar = 1.2 + 0.2. There is no secondary bar in 
this model. 

The amplitude spectra of the relative density perturba- 
tions (see e.g. Masset & Tagger, 1997; Rautiainen & Salo, 1999) 
(Fig-EJ show that the bar mode is not the only one in the disk, 
but there are also slower modes. The strongest of these modes, 
hereafter the S 1 mode, has an overlap of resonance radii with the 
bar: the corotation radius of the bar is approximately the same as 
the inner 4/ 1 -resonance radius of the slower mode (at T - 55 the 
RcR of the bar and the inner 4/1 resonance radius of the S 1 mode 
are both about 4.6 kpc). This resonance overlap does not seem 
to be a coincidence: when the amplitude spectra from different 
time intervals are compared, one can see that both the bar and 
the S 1 modes slow down so that the resonance overlap remains 
(see Fig. O. Furthermore, this reson ance overlap was t h e mos t 
common case in the simulations of iRautiainen & Salol (I1999I) . 
Also, the SI mode has a strong m - 1 signal and a maximum 
near its corotation at 7.1 kpc. The bar mode is also seen as a 
strong signal in the m - 4 spectrum, but only inside CR - the 
spiral part seems to be almost pure m = 2 mode. Altogether, the 
signals with m > 2 tend to be much weaker than features seen in 
m - I and m = 2 amplitude spectra. 

We have also tried to reconstruct the shapes of the modes 
seen in the amplitude spectra. This was done by averaging the 
surface density in coordinate frames rotating with the same an- 
gular velocities as the modes. No assumptions were made about 
the shapes of the modes. On the other hand, one should take 
these reconstructions with some caution, because the evolution 
of the two modes, the effect of slower (but weaker) modes, and 
short-lived waves may affect them. The results for the bar and the 
SI mode at the time interval T = 50-60 are shown in Fig.|4] The 
mode Qp = 47 km s^'kpc"' clearly shows the bar and symmetri- 
cal spiral structure that forms an Ri outer ring or pseudoring. By 
the T = 50-60 interval, the density amplitude of the bar mode is 
about 15-20 per cent in the outer ring region, where the maxima 
and minima have roughly the same strength. On the other hand, 
by r = 50 - 60, the mode Qp = 31 km s 'kpc ' is clearly lop- 
sided, which is not surprising considering the signal seen in the 
m - I amplitude spectrum. There is a minimum with an ampli- 
tude of about 30% and a maximum of about 15% aiR ~ 7 kpc, 
which corresponds to the CR of the SI mode. Earlier, at T ^ 20, 
the S 1 mode does not have the m - \ characteristic but exhibits a 
multiple-armed structure beyond its CR, accompanied by a clear 
signal in the m - 3 amplitude spectrum. 

3.3. The morphological changes in the gas component 

The amplitude spectra for the gas component at the interval 
r = 50 - 60 are also shown in Fig. |3] Due to fewer particles, 
they include more noise, but otherwise they are quite similar. In 



addition to the bar mode, the S 1 mode is also seen, but now it is 
more conspicuous in the m = 1 spectrum. 

The result of having several modes is the quite complicated 
evolution of the model (see Fig. |5]i: at different times, the mor- 
phology of the outer gaseous disk can be described as R]R'2, R'2, 
/?j or just as open spiral arms, which can sometimes be followed 
over 400 degrees. There is no evolutionary trend between the 
morphological stages, since they all appear several times during 
the model time span. The shape of the inner ring also changes by 
being sometimes more elongated or even consisting of tightly 
wound pair of spiral arms. On the broader sense, the overall 
Hubble stage of the model stays the same for several Gyrs. 

Although the slow modes in the stellar component can be 
clearly seen outside the bar radius (about 4 kpc), they become 
pronounced in the gas from R x 6 kpc. To study their effect on 
the gas morphology, we selected gas particles located at the an- 
nulus 7 < R < 10 kpc and calculated their number within every 
5°-sector along 6. Such density profiles were built for 301 mo- 
ments from the interval T=30-60 {T x 3-6 Gyr) with a step 
AT - 0.1 (~ 10 Myr). Earlier stages were not considered, be- 
cause then the pattern speed of the bar was changing so fast that 
it complicated the analysis. At every moment the distribution of 
gas density along was approximated by one-fold (m =1), two- 
fold (m = 2), and four-fold {m - 4) sinusoidal wave: 

cr = (To + A„ cos(m6 + 0,,,), (2) 

where cr is the gas density in a segment, ctq is the average density 
in the annulus, 0,„ and A„, are the phase and amplitude of the 
corresponding sinusoidal approximation, respectively. 

Figure |6] demonstrates the motion of maxima in the distri- 
bution of gas particles along 0. We made the density profiles in 
the reference frame co-rotating with the bar, whose major axis 
is always oriented in the direction 6-0°. Azimuthal angle 6 
is increasing in the sense of the galactic rotation, so the sup- 
posed position of the Sun is about - 315°. To illustrate the 
motion of density crests, we selected two intervals T=35.5-37.5 
and T=52.5-54.5 with a high amplitude of density perturbation. 
These density profiles indicate the motion of density maxima in 
the opposite direction to that of galactic rotation (i.e. they actu- 
ally rotate more slowly than the bar), which means an increase 
in the phase 0,„ of the sinusoidal wave (Eq.|2]i. 

Figure |7] exhibits the variations in the phase <p„ and ampli- 
tude A„, of the sinusoidal wave at the time intervals T=30-40, 
40-50, and 50-60. The subscripts 1 and 2 are related to the one- 
and two-fold sinusoids. Rotation of the density maxima causes 
the sharp changes in the phase when it achieves the value of 
<p - 360°, and at the new turn its value must fall to zero. These 
changes enable us to accurately calculate the mean values of the 
periods for the propagation of the sinusoidal waves, which ap- 
peal" to be f 1 = 3.3 + 0.4 and f 2 = 1 -5 ± 0.4. Remember that we 
study the density oscillations in the reference frame co-rotating 
with the bar, so the period P of beating oscillations between the 
bar and slow modes is determined by the relation: 



m(D.h - 

The periods. Pi and P2, appear to correspond to slow modes 
rotating with the pattern speeds Q = 28 + 2 km s ' kpc ' and 
Q = 26 ± 6 km s"' kpc respectively. It is more convenient 
to use simulation units here. The transformation coefficient be- 
tween them and (km s"' kpc"') is A: = 9.77, and the value of Q.h 
is Q/, = 4.8 s.u. The m- A wave manifested itself in two density 
maxima separated by the angle A0 » 90° (Fig. |6] right panel). 
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Fig. 4. The reconstructed modes in the stellar component (see text) for T = 50-60 time interval. The enhanced density compared to 
the azimuthally averaged profile at each radius is shown. The shades of gray (darker corresponds to higher surface density) have been 
chosen to emphasize the features. The circles in the bar mode indicate ILR (1.4 kpc), CR (4.6 kpc), and OLR (8.1 kpc), whereas the 
inner 4/1 (4.6 kpc) and CR (7.1 kpc) are shown for the mode SI. 



The analysis of phase motion of four-fold sinusoid reveals the 
period f 4 - 0.81 + 0.15, which corresponds to slow mode rotat- 
ing with the speed Q,/ = 28+4 km s"' kpc"' (Eq.[3]). Probably, 
it is mode Q = 28 + 4 km s"' kpc"' that causes the strong vari- 
ations in gas density with the periods Pi = 3.3, P2 - 15, and 
P4 - 0.8 when it works as m = 1, m = 2, and m - A density 
perturbations, respectively. This mode is well-defined in the gas 
and star power spectra made for the interval T=50-60 (Fig. [3]). 

Let us have a look at the amplitude variations (Fig. |7]i. The 
highest value of A2 equal A2 - 200 (particles per 5°-sector) is 
observed at the time T - 36.0 (left panel). On the other hand, Ai 
achieves its highest value of Ai = 220 at the time T - 56.5 (right 
panel). Amplitude A4 reaches its maximum value of A4 = 180 
at the time interval T = 53 - 55. Thus, the highest values of the 
amplitudes Ai, A2, and A4 are nearly the same. 

Figure |6] (left panel) indicates the growth of the amplitude 
of m=2 perturbation under a specific orientation of the density 
clumps. The amplitude of the sinusoidal wave is at its maximum 
at the moments T - 36.0 and 37.5 when the density clumps 
are located near the bar's minor axis, 9 - 90° and 270°. This 
growth is also seen in Fig. |7] (left panel) for the interval T=30- 
40: the amplitude A2 is at its maximum at the moments when 
(p2 ~ 180°. This phase corresponds to the location of maxima of 
OT = 2 sinusoid at 6 = 90° and 6 = 270° (Eq.|2]i. 

Our analysis revealed slight variations in the speed of the 
strongest slow mode, and they depend on its orientation with 
respect to the bar: Fig. Q (left panel) shows that the tilt of the 
phase curve, 02(f)> is variable. We can see that the slow mode 
rotates a bit faster when 02 ~ 180° (density clumps are near the 
bar's minor axis) and more slowly when <p2 - Q or 360° (the 
clumps are near the bar's major axis). Probably, the variations in 
the speed of the slow mode are connected with the change in the 
form of the density crests due to tidal interaction between the bar 
mode (barH-7?i-ring) and the slow mode. 



4. Kinematics of gas particles. Comparison with 
observations 

4. 1 . Momentary and average velocities 

We start our kinematical study with the interval T=50-60 (5- 
6 Gyr in physical time). At this period the bar rotates with a 
nearly constant pattern speed of Q.t, = 47 km s"' kpc"' which 
simplifies the analysis. The interval T=50-60 also provides the 
best agreement between the model and observed velocities. 

We determined the positions and velocities of gas particles 
at 101 moments separated by the step Ar = 0.1. For each mo- 
ment we selected gas particles located within the boundaries of 
the stellar gas complexes and calculated their mean velocities 
and velocity dispersions. To determine the positions of the com- 
plexes, we need to choose the position angle of the Sun with 
respect to the bar elongation, Ob- In this section we adopted the 
value of 6h = 45°, which gives the best fit between the model 
and observed velocities. 

Figure |8]shows the variations in the mean residual velocities, 
Vr and Vo, calculated for five complexes at different moments. 
The residual velocities were computed as differences between 
the model velocities and the velocities due to the rotation curve. 
It is clearly seen that the momentary velocities oscillate near the 
average values within the limits of ~ +20 km s"' . Two processes 
are probably responsible for these oscillations. The first is the 
slow modes that cause a quasi-periodic low in the velocity varia- 
tions. The second process is likely connected with the short-lived 
perturbations, e.g. from transient spiral waves in the stellar com- 
ponent. The averaging of velocities over long time interval re- 
duces the influence of slow modes and occasional perturbations. 

Table |2] represents the average values of the momentary 
residual velocities, Vr and Ve, calculated over 101 moments. It 
also gives the average values of velocity dispersions, and 0^, 
and the average number of particles « in the complexes. Since 
the bar has two tips, we calculated velocities for two opposite 
positional angles, 61, = 45° and Ot, - 225°, and used their mean 
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Fig. 5. The gas morphology at selected times. The bar is vertical in all frames, whose width is 20 kpc. 



Table 2. Model residual velocities averaged on interval T=50-60 



Region 


Vr 


CTr 


Vo 


0"9 


n 




km s"' 


km s"' 


km s"' 


km s"' 




Sagittarius 


8.5 


7.2 


0.1 


5.9 


70 


Carina 


7.5 


7.6 


-2.0 


6.6 


158 


Cygnus 


6.8 


10.1 


8.2 


6.5 


108 


Local System 


6.8 


11.7 


4.8 


6.7 


112 


Perseus 


-12.5 


11.9 


-2.9 


6.5 


70 



values. The averaged residual velocities are determined with the 
errors of 0.4-1.4 km s '. The relatively low level of the errors is 
due to the large number of moments considered (N=101). 

When comparing Tables |2]and [1] one can see that our model 
reproduces the directions of the radial and azimuthal compo- 
nents of the residual velocities in the Perseus and Sagittarius 
regions and in the Local System. We succeed in the Sagittarius 
region where our model reproduces the observed velocities with 
the accuracy 1.4 km s"'. Unfortunately, in the Perseus region 
the model residual velocity \Vr\ is too high, and the difference 
between the model and observed velocities achieves 5.8 km s"' 



there. Our model can also reproduce the positive Vr velocity in 
the Local System, which deviates only 1.5 km s ' from the ob- 
served one. 

We now consider the mean difference between the model and 
observed velocities AV calculated for the radial and azimuthal 
components: 

.^2 _ i:i{(^-^«obs)' + c^-^«obs)1 ,,, 

where A; is a number of complexes. The value of AV computed 
for three complexes (the Sagittarius and Perseus regions and the 
Local System) equals AV - 3.3 km s"'. Another situation is ob- 
served in the Carina and Cygnus regions where we cannot even 
reproduce the direction of the observed residual velocities. The 
mean value of the velocity deviations achieves AV - 13.3 km 
s"' there. 

To demonstrate the distribution of the average velocities on 
the galactic plane, we divided the area (-10 < x < +10, 
-10 < y < +10 kpc) into small squares of a size 0.250 x 0.250 
kpc. For each square we calculated the average values of the 
residual velocities of gas particles. Then we averaged residual 
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Fig. 7. Variations in the phase (black curve) and amplitude A (gray curve) of the sinusoids that approximate the distribution of 
gas particles located at the distances 7 < /? < 10 kpc along 0. Subscripts 1 and 2 are related to the one- and two-fold sinusoids, 
respectively 
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Fig. 8. Variations in the mean velocities of gas particles located within the boundaries of the stellar-gas complexes. The left panel is 
related to the radial component Vr and the right one to the azimuthal one Vg. 



velocities over 101 moments for the interval T=50-60. The av- 
erage residual velocities in squares are shown in Fig.|9] We de- 
picted only squares that contain high enough number of parti- 
cles, n > m/2, where « is the number of particles accumulated in 
a square over 101 moments but 7n is their number averaged over 
all squares, m - 463. 

In Paper I we have built similar figures for models with an- 
alytical bars. Two different moments were considered: when the 
broken rings (pseudorings) were observed and when they trans- 
formed into pure rings. The pseudorings and pure rings created 
different kinematical pictures. We connected the main kinemati- 
cal features of the pseudorings with the gas outflow and those of 
the pure rings with the resonance. The distribution of the neg- 
ative and positive velocities obtained for N-body simulations 
(Fig. |9]l strongly resembles that of the pseudorings in models 
with analytical bars, giving support to the "averaging process" 
adopted here. This similarity suggests there is gas outflow in the 
present model (see also Sect. 4.4). 

4.2. Velocities in the complexes under different values of the 
solar position angle Ob 

We studied the dependence of the average residual velocities 
Vr and Ve on the solar position angle Of,. Figure [TOhb shows 



5 curves that demonstrate the velocity changes in 5 complexes. 
The sharpest changes in the radial velocity Vr are observed in 
the Local System and in the Cygnus region, and the radial veloc- 
ities in the other complexes depend only weakly on the choice 
of Of,. As for the azimuthal component, the strongest changes 
can be seen in the Sagittarius, Carina, and Perseus regions, but 
the velocity changes are modest in other complexes. Practically 
speaking, the optimal value of 61, provided the best agreement 
between the model and observed velocities is determined by the 
radial velocity in the Local System and by the azimuthal velocity 
in the Sagittarius region. These velocities achieve their observed 
values of Vr = 5.3 and Vg - -1 km s"' under 61, - 43° and 
6b - 48°, respectively. 



We now consider the sum of square differences between 
the model and observed velocities, obtained for the radial 
and azimuthal components under different values of 9b. Figure 
[TOb shows the ;t'^-function computed for three complexes: the 
Perseus and Sagittarius regions and the Local System. It is 
clearly seen that achieves its minimum values at the inter- 
val 9b > 40°. We chose 9b = 45° as the optimal value because 
it reproduces the observational velocity Vg = -1 km s ' well 
in the Sagittarius region. Models with analytical bars in Paper I 
gave the same result. 
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Fig. 9. Distribution of the negative and positive average residual velocities calculated in squares. The squares with positive velocities 
are shown in black, while those with negative ones are given in gray. Only squares that satisfy the condition n > m/2 are depicted. 
The left panel represents the radial velocities, while the right one shows the azimuthal ones. It also demonstrates the boundaries of 
the stellar-gas complexes. The position angle of the Sun is supposed to be 9b = 45°. The bar is oriented along the Y-axis, the galaxy 
rotates clockwise, and a division on the X- and F-axis corresponds to 1 kpc. 



4.3. Analysis of periodicity in oscillations of the momentary 
velocities 

Now we approximate the oscillations in the radial and azimuthal 
components of the momentary velocities, Vr and Ve (Fig. [8]), by 
the sinusoidal law: 



Vr(oi Ve) = Ai sm(2nT/P) + A2 cos(2nT/P), 



(5) 



where P is a period of oscillations, Aq = ^A^ + is an ampli- 
tude of oscillations, and T is time counted from To - 50. 

We use the standard least square method to solve the system 
of 101 equations, which are linear in the parameters Ai and A2 
for each value of nonlinear parameter P. We then determine the 
value of P that minimizes the sum of squared normalized resid- 
ual velocities x^- Figure [TT] presents the -curves built for the 
oscillations of the radial velocity in 5 complexes, but the curves 
made for the azimuthal velocities have no conspicuous minima. 
It is clearly seen thatjt'^-curves demonstrate deep minima in the 
Cygnus and Perseus regions and in the Local System. These 
minima correspond to the best periods in approximating the ve- 
locity oscillations that have the following values: P - 2.7 + 0.4 
in the Cygnus region, P = 2.9 + 1.0 in the the Local System, and 
P = 1 .6+0.2 in the Perseus region. We have already obtained pe- 
riod P = 1 .5 when studying density oscillations on the galactic 
periphery (Sect. 3.3). Probably, the strongest slow mode Q = 28 
km s~' kpc"' is also responsible for the velocity oscillations: the 
beating oscillations between the bar mode and a two-armed pat- 
tern rotating with the speed Q = 28 km s~' kpc ' must have 
the period of f = 1.6 and those calculated for one-armed per- 
turbation have a period of P = 3.2. (Eq. [3]l. Some of the small 
differences between the pattern speeds derived from the ampli- 
tude spectra and those obtained from kinematical analysis may 
be due to tidal interaction in the stellar component between the 
bar and slow modes. 



4.4. Evolutional aspects of kinematics at the time interval 
7=30-60 

Let us compare the average residual velocities calculated for 
different time intervals T=30-40, 40-50, and 50-60 (Tables |2] 



[3] and |4|. Generally, most changes in the residual velocities 
do not exceed 4.0 km s ' and are likely caused by occasional 
perturbations. On the other hand, radial velocities Vr in the 
Local System and in the Cygnus region demonstrate the ongoing 
growth, which can be connected with the evolution of the outer 
rings. 

Figure [12] shows the surface density of gas particles aver- 
aged in squares at different time intervals. The average density 
was calculated in the reference frame that rotates with the speed 
of the bar. The light-gray, dark-gray, and black colors represent 
squares containing the increasing number of particles, n > m/2, 
n > m, and n > 2m, respectively, where n is the number of par- 
ticles accumulated in a square over 101 moments and m is their 
number averaged over all squares, m = 463. It is clearly seen 
that the major axis of the outer ring R2 changes its orientation: 
it goes a ~ 20° ahead of the bar at the interval T=30^0, but 
this angle increases to or ~ 45° at the intervals T=40-50 and 
T=50-60. Moreover, the outer ring changes its morphology: we 
can identify two outer rings of classes R[ and R2 at the interval 
T=30-40, while there is only one outer ring with an intermediate 
orientation of a w 45° at the intervals T=40-50 and 50-60. Its 
shape becomes rounder at the interval T=50-60. 

Let us consider more thoroughly the distribution of gas par- 
ticles at the interval T=50-60 (Fig.[T2li. It is clearly seen that the 
surface density of gas particles at the distance range of /? = 6 - 9 
kpc is nearly twice the average density all over the disk m. The 
density perturbation inside the outer ring can be approximated 
by two spiral arms with a pitch angle of / = 6 ± 1 ° . The density 
perturbation inside them reaches to 100 per cent with respect to 
the average gas density in the disk. This is considerably larger 
than the density perturbation seen in the stellar component (15- 
20 per cent). 

Figure [13] shows the profiles of the surface density of gas 
particles averaged at the different time intervals. We can see the 
growth of the density hump at the distance of R x 7 kpc, which 
indicates the growth of the outer ring. In contrast, the hump at 
R X 3 kpc is decreasing, which reflects the weakening of the 
inner ring. At the interval T - 50-60, the maximum in the gas 
density distribution is located at the distance R = 1.3 kpc, which 
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Fig. 10. Dependence of the average residual velocities, Vr (a) 
and Vg (b), and the -function (c) on the solar position angle Oi,. 
The curves calculated for different complexes are shown by dif- 
ferent lines. The -function was computed for three complexes: 
the Perseus and Sagittarius regions and the Local System. 

is just in the middle between the outer 4/1 resonance (6.4 kpc) 
and the OLR (8. 1 kpc) of the bar. 

Tables |2l [3] H also represent the velocity dispersions of gas 
particles in the stellar gas complexes. We can see that their av- 
erage values stay at nearly the same level ofo^ = 9.7 + 0.1 and 

- 6.3 + 0.2 km s ' during the period T=30-60. The maxi- 
mum growth, which does not exceed ~ 20 per cent, is observed 
in the Perseus region. The model velocity dispersions somewhat 
exceed the observed values derived for OB-associations in the 
stellar-gas complexes, cr^ gbs ~ ^-^ "^eobs ~ ^--^ ^ 
but this difference is below 30 per cent. 

5. Conclusions 

We have presented N-body simulations that reproduce the kine- 
matics of OB-associations in the Perseus and Sagittarius regions 
and in the Local System. The velocities of gas particles averaged 
over large time intervals (1 Gyr or 8 bar rotation periods) repro- 
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Fig. ll.;i'^-functions built for studying periodicity in the oscilla- 
tions of the radial velocities, Vr, in 5 complexes. The minima on 
the curves must correspond to the best periods in approximation 
of the velocity oscillations. 

Table 3. Model residual velocities averaged on interval T=30-40 



Region 




O-R 


Vo 




n 




km s"' 


km s"' 


km s"' 


km s"' 




Sagittarius 


10.2 


8.4 


1.8 


6.9 


84 


Carina 


10.4 


8.7 


1.2 


7.8 


184 


Cygnus 


-0.7 


9.9 


10.8 


5.4 


87 


Local System 


-1.0 


11.6 


7.7 


5.9 


90 


Perseus 


-11.3 


9.7 


-2.5 


5.6 


83 


ible 4. Model residual velocities averaged on interval T=40 


Region 


Vr 


o-R 


Ve 




n 




km s"' 


km s"' 


km s"' 


km s"' 




Sagittarius 


8.9 


7.4 


0.8 


5.7 


69 


Carina 


9.1 


7.4 


-1.0 


6.8 


164 


Cygnus 


3.9 


104 


11.4 


5.7 


126 


Local System 


5.6 


12.5 


8.3 


6.7 


121 


Perseus 


-15.2 


11.1 


-1.9 


5.5 


50 
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Fig. 12. The surface density of gas particles averaged in squares at the time intervals T=30^0, 40-50, and 50-60. The light- 
gray, dark-gray, and black colors represent squares containing the increasing number of particles: n > m/2, n > m, and n > 2m, 
respectively. The position of the Sun is shown by the specific symbol. The bar is oriented along the Y-axis, the galaxy rotates 
clockwise, and a division on the X- and Y-axis corresponds to 1 kpc. 




2 4 6 8 10 R, kpc 



Fig. 13. Profiles of the surface density of gas particles averaged 
at the intervals T=30-40, 40-50, and 50-60. 

duce the directions of the observed velocities in these regions. 
The mean difference between the model and observed velocities 
calculated for the radial and azimuthal components is AV = 3.3 
km s ' there. 

The galactic disk in our model includes two subsystems. 
The behavior of the stellar subsystem is modeled by 8 mil- 
lion gravitating collisionless particles. The stellar disk quickly 
forms a bar. Its original pattern speed is quite high, but it 
first quickly decreases and then moves to a slow decrease with 
n «i 50 km s 'kpc ' for several Gyrs. With our favored value 
of the solar distance, Rq = 7.5 kpc, this sets us close to the 
OLR (RoLR - 8.1 kpc). This agrees with studies of local stellar 
velocity distribution (Dehnen, 2000; Fux, 2001; Minchevet al., 
l2009l) . although they tend to set the OLR slightly inside Rq. 
The optimal value of the solar position angle Oh providing the 
best agreement between the model and observed velocities is 
0b = 45 + 5°. The bar is quite long (Rbar ~ 4.0 kpc), but 
both its size and orientation are consis tent with the parame- 
ters derived from infrared o bservations dBeniamin at al.L 120051: 
ICabrera-Lavers et al.Ll2007l) . 



The stellar disk also creates an outer ring of class Ri rotating 
with the pattern speed of the bar, and the corresponding density 
perturbation amounts to 15-20 per cent of the average density at 
the same distance. Besides the bar, the stellar disk includes sev- 
eral slow modes. The strongest of these rotates with the pattern 
speed of Q w 30 km s 'kpc ' and is often clearly lopsided. 

The gas subsystem is modeled by 40 000 massless particles 
that move in the potential created by the stellar particles (and 
analytical bulge and halo) and can collide with each other in- 
elastically. The gas disk forms an outer ring that exhibits quasi- 
periodic changes in its morphology because it has several modes. 
One can identify elements of Ri- and /?2-morphology, and the 
outer ring can often be classified as /^i/Jj. The gas density per- 
turbation inside the ring can be approximated by two spiral arms 
with the pitch angle of / = 6 + 1°. 

The models with analytical bars (Paper I) reproduced the 
residual velocities well in the Perseus and Sagittarius regions. 
We explained this success by the resonance between the rela- 
tive orbital rotation of the bar and the epicyclic motion. The 
Sagittarius region must be located slightly inside the OLR where 
resonance orbits are elongated perpendicular to the bar, whereas 
the Perseus region must lie outside the OLR where periodic or- 
bits are oriented along the bar However, models with the ana- 
lytical bar failed dramatically with the Local System where they 
yielded only negative radial velocities -24 < Vr < -16 km s"', 
whereas the observed value is +5.3 km s"'. The success of N- 
body simulations with the Local System is likely due to the 
gravity of the stellar /?i-ring, which is omitted in models with 
analytical bars. 

To study the effects of the gravity of the 7? i -ring we con- 
struct more simple models with a "time averaged bar potential". 
This was done by calculating the average density distribution in 
the frame rotating with the bar This process should average out 
most of the effect of slower modes and leave bar and the /?i-ring 
that corotates with the bar The preliminary study shows that mo- 
mentary velocities in these models are in a good agreement with 
the average velocities in the present N-body simulation. The de- 
tailed description of these models will be done in our next paper 

To simplify the analysis at this point we are forced to ignore 
a lot of processes which are important at such long time inter- 
val as 6 Gyr We do not consider the accumulation of gas at the 
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galactic center, the transitions between the gas and stellar sub- 
systems, resonant interaction between the bar and halo, or the 
minor mergers and satellite accretion. Considering the effects of 
these processes may be done in a later phase. 
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